Homework 4 , Spring 2022
Name: Anitha Ganapathy
Email: aganapa@iu.edu
Organizing the imports
In [1]:
# Import the necessary libraries
from PIL import Image
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# import matplotlib.mlab as mlab
import librosa
import soundfile as sf
from IPython.display import Audio
import cmath, math
import scipy.io # to read .mat files
import seaborn as sns
from numpy.core.fromnumeric import size
import scipy
from google.colab import drive
drive.mount('/content/drive')
In [2]:
dirpath = '/content/drive/MyDrive/MLSP_data/'
!ls $dirpath
P1: KNN Source Separation¶
In [3]:
trs_file = dirpath + 'trs.wav'
trn_file = dirpath + 'trn.wav'
trs, sr = librosa.load(trs_file , sr = None)
trn, sr = librosa.load(trn_file , sr = None)
In [4]:
trs.shape, trn.shape
Out[4]:
Original S and N i.e piano and ocean signals¶
In [5]:
Audio(trs, rate = sr)
Out[5]:
In [6]:
Audio(trn, rate = sr)
Out[6]:
In [7]:
hop_length = 512
In [8]:
S = librosa.stft(trs, n_fft=1024, hop_length = hop_length, window='hann')
N = librosa.stft(trn, n_fft=1024, hop_length = hop_length, window='hann')
In [9]:
S.shape , N.shape
Out[9]:
Matrix G, the magnitude spectrogram of trs.wav+trn.wav.¶
In [10]:
S_abs = np.abs(S)
N_abs = np.abs(N)
In [11]:
G = S_abs + N_abs
G.shape
Out[11]:
In [12]:
Audio(librosa.istft(G), rate = sr)
Out[12]:
Ideal Binary Masks (IBM),¶
B
$$ \boldsymbol{B}_{f, t}= \begin{cases}1 & \text { if } \boldsymbol{S}_{f, t} \geq \boldsymbol{N}_{f, t} \\ 0 & \text { otherwise }\end{cases} $$The hope is that
S ≈ B ⊙ G
In [13]:
B = (S_abs > N_abs).astype(int)
B.shape
Out[13]:
In [14]:
B
Out[14]:
Load x_nmf.wav and convert it into a complex-valued spectrogram X and its magnitudes Y¶
In [15]:
x_nmf_file = dirpath + 'x_nmf.wav'
x_nmf, sr = librosa.load(x_nmf_file , sr = None)
In [16]:
x_nmf.shape
Out[16]:
In [17]:
Audio(x_nmf, rate = sr)
Out[17]:
In [18]:
X = librosa.stft(x_nmf, n_fft=1024, hop_length = hop_length, window='hann')
In [19]:
X.shape
Out[19]:
In [20]:
plt.specgram(x_nmf ,Fs = sr)
plt.suptitle("x_nmf wav, signal in frequency domain using the spectrogram ", y=1.05, fontsize=18)
plt.xlabel('Time')
plt.ylabel('Frequency')
plt.show()
In [21]:
Y = np.abs(X)
Y.shape
Out[21]:
In [22]:
G.shape
Out[22]:
In [23]:
def knn_source_separation(Y, G, k,B, X):
print("Y shape :", Y.shape) # Y shape : (513, 131)
print("G shape :", G.shape) # G shape : (513, 987)
euc_dist = np.zeros((Y.shape[1], G.shape[1]))
# we compare each column of Y with all columns in G to find out knearest neighbors.
for i in range(Y.shape[1]): # (columns : 131)
for j in range(G.shape[1]): # (columns : 989)
euc_dist[i][j] = np.sqrt(np.sum((Y[:, i] - G[:, j]) ** 2))
sorted_distances = euc_dist.argsort()
knn_indexes = sorted_distances[:,0:k] # k indexes with shortest distance
print("knn_indexes shape : ", knn_indexes.shape) # knn_indexes shape : (131, k)
# print(knn_indexes)
# D = IBM of the test signal
D = np.zeros((Y.shape[1],Y.shape[0])) # Y shape : (513, 131)
# B shape : (513, 989)
# By using the nearest neighbor indices, I can collect the IBM vectors
# from B, i.e. {B:,i1 ,B:,i2 , · · · ,B:,ik}.
for m in range(knn_indexes.shape[0]): # 131 for every row
b_val = np.zeros((513,1))
for n in range(knn_indexes.shape[1]): # 20 for every column index
# b_val.append(np.median(np.sort(B[:, knn_indexes[m][n]])))
col_index = int(knn_indexes[m][n])
median_value = B[:,col_index]
median_value = median_value.reshape((median_value.shape[0],1))
b_val = np.hstack((b_val,median_value))
median_k = np.median(b_val[:,1:], axis = 1)
D[m,:] = median_k
print(f'D.T shape : {D.T.shape}')
print(f'X shape : {X.shape}')
S_recov = np.multiply(D.T,X)
print(f'S_recov shape : {S_recov.shape}')
return S_recov
In [24]:
source_recov = knn_source_separation(Y, G, 12,B, X)
X_nmf recovered source signal¶
In [25]:
source_actual = librosa.istft(source_recov, hop_length = hop_length)
Audio(source_actual, rate = sr)
Out[25]:
In [25]:
P2: Motor Imagery¶
In [26]:
!ls $dirpath
In [27]:
eeg_mat = scipy.io.loadmat(dirpath+'eeg.mat')
In [28]:
print(eeg_mat.keys())
# print(mat.get('L'))
In [29]:
train_x = eeg_mat.get('x_train')
train_y = eeg_mat.get('y_train')
test_x = eeg_mat.get('x_te')
test_y = eeg_mat.get('y_te')
In [30]:
train_x.shape, train_y.shape, test_x.shape, test_y.shape
Out[30]:
In [31]:
print(train_x[:,0,:].reshape(-1).shape)
train_x[:, 0, :].shape
Out[31]:
Blackman window¶
In [32]:
def blackman_window(data, N, hop_length):
blackman = np.blackman(N)
new_data = []
hop_counter = 0
while hop_counter + N <= len(data):
frame = data[hop_counter: hop_counter+N]
vector = np.multiply(frame.T, blackman)
new_data.append(vector)
hop_counter = hop_counter + hop_length
return np.array(new_data).T
In [33]:
def stft_for_mu_freq(data):
stft_mu_vector = []
dft_eeg_sig = np.fft.fft(np.eye(64))
for i in range(3):
channel_window = blackman_window(data[:,i], 64, 48)
stft_mu = np.abs(np.dot(dft_eeg_sig, channel_window))[2:7,:].reshape(-1,1)
stft_mu_vector.append(stft_mu)
return np.array(stft_mu_vector).reshape(-1,1)
Creating the Z matrix
In [34]:
Z_train = []
Z_test = []
for i in range(train_x.shape[2]):
Z_train.append(stft_for_mu_freq(train_x[:,:,i]))
for i in range(test_x.shape[2]):
Z_test.append(stft_for_mu_freq(test_x[:,:,i]))
print(len(Z_train), len(Z_test))
print(len(Z_train[0]), len(Z_test[0]))
In [35]:
Z_train = np.array(Z_train).T.reshape(225,112)
Z_test = np.array(Z_test).T.reshape(225,28)
print(Z_train.shape, Z_test.shape)
In [36]:
print(Z_train.shape, Z_test.shape)
In [37]:
def matrix_A_calc(l,M):
np.random.seed(3232)
A = np.random.uniform(-1,1,l)
for i in range(len(M) - 1):
A = np.vstack((A,np.random.uniform(-1,1,l)))
sum_A = [sum(A[i]) for i in range(len(M))]
sum_A = np.array(sum_A)
inv_sum_A = 1/sum_A
B = A * inv_sum_A[:, np.newaxis]
return B.T
def matrix_Y_calc(A,Z):
Y = np.dot(A,Z)
Y_sign = np.sign(Y)
return Y,Y_sign
def hamming_dist(s1,s2):
return np.count_nonzero(s1 != s2)
def distance(Y,Y_test):
distance = np.zeros((28,112))
for i in range(distance.shape[0]):
for j in range(distance.shape[1]):
distance[i][j] = hamming_dist(Y_test[:,i],Y[:,j])
sorted_distance = distance.argsort()
final_index = np.zeros((Y_test.shape[1],Y.shape[1]))
for i in range(Y_test.shape[1]):
for j in range(Y.shape[1]):
index = sorted_distance[i][j]
final_index[i][j] = train_y[index,0]
return final_index
KNN algo for the given k = 12¶
In [38]:
l_list =[]
acc_list = []
k = 12
for i in range(2,21,2):
l = i
A = matrix_A_calc(l,Z_train)
Y,Y_sign = matrix_Y_calc(A,Z_train)
Y_test,Y_test_sign = matrix_Y_calc(A,Z_test)
knn_dist_sorted = distance(Y_sign,Y_test_sign)
kNearestNeighbours = knn_dist_sorted[:,0:k]
final_y_test = np.zeros((test_y.shape[0],1))
for p in range(0,28):
final_y_test[p] = np.median(kNearestNeighbours[p,:])
count = 0
for p in range(0,28):
if final_y_test[p] == test_y[p]:
count+=1
accuracy = count/28
l_list.append(i)
acc_list.append(accuracy)
LSH and accuracy values¶
In [39]:
df = pd.DataFrame(list(zip(l_list, acc_list)), columns= ['LSH', 'Acc'])
df
Out[39]:
In [40]:
# plot the bar graph
fig = plt.figure(figsize=(8,5))
ax = fig.add_axes([0,0,1,1])
plt.bar(l_list,acc_list, color ='lightblue',width = 0.9)
plt.plot(l_list,acc_list, color ='green')
plt.title("Bar plot of lsh vs the accuracy", fontsize = 20)
plt.suptitle("k value is : " +str(k), fontsize = 15)
plt.xticks(l_list)
plt.xlabel("LSH values", fontsize = 18)
plt.ylabel("accuracy", fontsize = 18)
plt.show()
P3: Multidimensional Scaling¶
In [41]:
!ls $dirpath
In [42]:
mat = scipy.io.loadmat(dirpath+'MDS_pdist.mat')
In [43]:
print(mat.keys())
print(mat.get('L'))
In [44]:
mds = np.array(mat.get('L'))
print("MDS Data shape: ",mds.shape)
In [45]:
mds[:5, :5]
Out[45]:
MultiDimensional Scaling (MDS) - M08, S05¶
The recipe for MDS
- Convert your pairwise distance matrix $M$ into a positive definite matrix $W$
- Do eigendecomposition on the matrix $W$
- First few vectors of $UV^{1/2}$ holds recovers the locations (coordinates)
In [46]:
mds.shape
Out[46]:
In [47]:
np.max(mds) , np.min(mds)
Out[47]:
In [48]:
plt.figure(figsize=(5,5))
c= plt.pcolormesh(mds[:15, :15])
plt.show()
plt.figure(figsize=(15,8))
c= plt.pcolormesh(mds[:, :])
plt.colorbar(c)
plt.title('color mesh of MDS data', fontweight ="bold")
plt.show()
In [49]:
np.mean(mds, axis=1).shape # row wise mean
Out[49]:
In [50]:
m_bar = np.mean(mds, axis=1)
m_bar.shape
Out[50]:
In [51]:
m_tilde = mds - m_bar
m_tilde.shape
Out[51]:
In [52]:
X_mat = m_tilde - np.mean(m_tilde, axis=0) # substraction of row wise mean
In [53]:
X_mat.shape
Out[53]:
In [54]:
def eigen_decomposition(X, num_eigenvectors):
max_epochs = 50
eigen_val = []
eigen_vec = []
for i in range(num_eigenvectors):
w = np.random.rand(X.shape[0], 1)
for epoch in range(max_epochs):
w = np.dot(X, w)
w = w * (1/np.sqrt(np.sum(w ** 2)))
s1 = np.sqrt(np.sum(np.dot(w.T, X) ** 2))
u1 = np.dot(X.T, w) /s1
X = X - (s1 * np.dot(w, u1.T))
eigen_val.append(s1)
eigen_vec.append(w.reshape(-1))
eigen_vec = np.array(eigen_vec)
return eigen_vec.T, eigen_val
In [55]:
eigenvectors, eigenvalues = eigen_decomposition(X_mat, 2)
In [56]:
print(eigenvectors)
print(eigenvalues)
In [57]:
np.diag(eigenvalues)
Out[57]:
In [58]:
matrix = np.dot(eigenvectors, np.sqrt(np.diag(eigenvalues)))
In [59]:
print(matrix.shape)
matrix
Out[59]:
In [60]:
plt.figure(figsize=(20, 4))
# c= plt.pcolormesh(mds[:15, :15])
c= plt.pcolormesh(matrix.T)
plt.colorbar(c)
plt.title('color mesh of MDS data', fontweight ="bold")
plt.show()
Report the map you recovered in the form of a scatter plot.¶
In [61]:
plt.scatter(matrix[:,0], matrix[:,1], )
plt.title('MDS map', fontsize =20)
plt.show()
plt.scatter(matrix[:,1], matrix[:,0], )
plt.title('Rotated map', fontsize =20)
plt.show()
Final: Convert the notebook to HTML¶
In [62]:
# %%shell
# jupyter nbconvert --to html /content/AG_MLSP_SP22_HW_4.ipynb